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We review studies of vibrational energy transfer in a molecular junction geometry, consisting of a 
molecule bridging two heat reservoirs, solids or large chemical compounds. This setup is of interest 
for applications in molecular electronics, thermoelectrics, and nanophononics, and for addressing 
basic questions in the theory of classical and quantum transport. Calculations show that system 
size, disorder, structure, dimensionality, internal anharmonicities, contact interaction, and quan¬ 
tum coherent effects, are factors that interplay to determine the predominant mechanism (ballis¬ 
tic/diffusive), effectiveness (poor/good) and functionality (linear/nonlinear) of thermal conduction 
at the nanoscale. We review recent experiments and relevant calculations of quantum heat transfer 
in molecular junctions. We recount the Landauer approach, appropriate for the study of elastic 
(harmonic) phononic transport, and outline techniques which incorporate molecular anharmonici¬ 
ties. Theoretical methods are described along with examples illustrating the challenge of reaching 
control over vibrational heat conduction in molecules. 
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I. INTRODUCTION 


Achieving control over energy transfer at the molecular scale is essential for the realization of molecular-based tech¬ 
nologies. Good thermal conductors are desired in electronic applications where it is necessary to efficiently remove 
excess heat generated by electronic heat dissipation m- On the other hand, materials which are poor conductors 
of heat are required for the design of efficient thermoelectric devices HIS]. Intramolecular vibrational energy redis¬ 
tribution (IVR) and intermolecular vibrational energy transfer are central processes in chemistry, fundamental to our 
understanding of chemical dynamics, specifically, reaction rates PH]. Moreover, in biomolecules, energy transfer 
processes determine stability, function, and regulation |5]. 

In this review, motivated by progress in probing heat transfer at the nanoscale, particularly, in self-assembled 
monolayers (SAMs) and nanostructured materials [9], we focus on the problem of vibrational heat flow through 
molecular junctions. After the presentation of representative experimental results, we describe quantum mechanical- 
based methodologies for the simulation of vibrational energy flow, focusing on the behavior of a single molecule. In 
the generic setup of interest a molecule is placed between two macroscopic contacts, labeled by L and R, see Fig. 
[^a) for a schematic illustration. These thermal baths are characterized by well-defined thermodynamic properties, 
essentially, their temperatures Tl and Tr, potentially creating large biases, (Tr — Tii)/{TL + Tr) ~ I. Considering 
transport of vibrational energy in the junction, basic quantities of interest are the heat current jq and the thermal 
conductance k = jq/AT (units W/K) with AT = Tr — Tr as the temperature difference between the bulk objects. 
Complementary to solid/molecule/solid junctions, heat transfer in molecular chains can by studied in a “bridged” 
configuration, by connecting a molecule at its two ends to relatively large side groups which act (approximately) as 
thermal baths, see Fig. [^b) for a scheme following Ref. m- In such experiments, a short laser pulse is applied to 
heat one end of the extended molecule. The heat transfer rate is then identified from optical-spectroscopy means, by 
measuring either the rate of energy loss from the excited unit, or energy gain at the other end. 

Besides its relevance to technologies such as molecular electronics [mill], phononics [131 m], thermoelectrics |S], 
spin caloritronics HU, N anoelectromechanics m, quantum optomechanics M and quantum information processing 
with hybrid quantum circuits |18j . the generic model displayed in Fig. [T a) is of basic interest for probing open 
quantum systems phenomena. Energy transfer processes in molecular junctions involve many-body interactions and 
quantum effects in a non-equilibrium situation, thus the setup can be employed to address fundamental questions in 
statistical mechanics, condensed matter physics, thermodynamics, and classical and quantum mechanics. Is energy 
transport in a given system ballistic or diffusive? What are the roles of many-body interactions in determining 
transport characteristics, especially when the system is driven beyond linear response? What are the signatures of 
quantum effects in molecular-scale heat transport? Can we design or predict a certain functionality from the model 
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FIG. 1: (a) Scheme of a molecular junction consisting of two solids (possibly made of different materials, gold and 
quartz) and a molecule of several units, typically functionalized at its ends to improve chemical bonding with the 
solids, e.g., alkane-thiols of N units. Heat current is flowing from the hot reservoir (e.g., the left one) to the 
other-cold side, (b) Vibrational energy flow in molecules can be measured using pump-probe laser spectroscopy tools, 
in liquid solutions at ambient temperatures. In this illustration, the molecule of interest (alkane chain) is connected 
at its two ends to large molecular groups acting as local heat reservoirs, for example, azulene and anthracene m- 
The azulene group is excited optically. The time for intramolecular vibrational energy transfer, through the bridge 
to the anthracene unit, is measured for different chain lengths, to identify the dominant transport mechanism. 


Hamiltonian? Other questions (not reviewed here), which have attracted significant attention, concern the emergence 
of the phenomenological-macroscopic Fourier’s law of heat conduction from first principles (classical and quantum) 
m, and the derivation of sufficient or necessary conditions for realizing negative differential thermal conductance 
and thermal rectification Hamlin]- The latter effect, also referred to as the thermal diode effect, features junctions 
in which the magnitude of the thermal current depends on the direction (polarity) of the temperature difference. 

In this review our main focus are theoretical frameworks for the calculation of vibrational energy flow in molecular 
junctions, and we contain ourselves to quantum mechanical treatments. Why concern with quantum effects? Molecular 
vibrational frequencies typically extend room temperature, huj > ksT, thus the quantum statistics is important; in 
classical simulations the average energy per phonon mode is ksT, with ks as the Boltzmann constant, resulting in 
an overestimate of e.g., the specific heat (thus the thermal conductivity), since < ksT. 

Because of the wide scope of the field of heat transfer, we leave out many interesting topics which were covered in 
recent reviews. Concepts, experiments, and computational approaches in nanoscale thermal transport were recently 
organized into several comprehensive works [smumn], with a focus on low-dimensional nanostructured materials 
such as graphene, carbon nanotubes, and Si nanowires, rather than molecular-level transport. Other relevant reviews 
cover control over energy dissipation at the nanoscale [T] , advances in thermoelectrics of semiconductor nanostructured 
materials [3, studies over the nature of vibrational energy flow in proteins [5|, and heating, heat transfer, and 
thermoelectricity at the nanoscale [23] . Concerning methodologies for calculating nanoscale phonon transport, recent 
reviews detailed classical simulations [lU EH ES] and quantum-mechanical methodologies m, specifically Green’s 
function-based approaches EHEHj- Other studies featured phononic devices mmi and phonon-assisted effects in 
molecular electronic conduction |29] . Obviously, the interaction of electrons with phonons is intrinsic for both charge 
and energy transfer processes in molecular systems m- Here we simplify the problem considerably by studying only 
the dynamics of vibrational-nuclear degrees of freedom. Even under this separation, the situation depicted in Fig. [l]is 
deeply complex given the interplay of quantum effects, many-body (phonon-phonon) interactions, far-from-equilibrium 
driving, and the diversity in molecular structures and contact geometries. 

Before discussing experimental results and theoretical methodologies, it is useful to recall two central-limiting energy 
transfer mechanisms which are well defined in macro-to-mesoscale systems: ballistic and diffusive. Ballistic transport 
corresponds to direct point-to-point propagation of energy, obeying the scaling I ex t, with t as the traveling time to 
cross the distance 1. It shows up in ordered structures in which the mean free path, the average distance traveled by 
phonons between successive scattering events, is larger than the size of the conductor. This is the case e.g., in some 
(up to ^m-length) nanotubes [3], and more relevant to our review, alkane chains of 5-25 units, at room temperature 
[ini[3lH36]. Diffusive energy transport is a multiple-step process, with energy “hopping” between sites, obeying the 
scaling I ex \/i. It is observed in bulk systems and disordered structures, e.g., peptides [SI [37], glasses, amorphous and 
doped nanostructures [3]. 

The distinction between ballistic and diffusive dynamics, and the notion of a mean free path, are natural to macro¬ 
scale and nanostructured systems, but in relatively short molecules as considered below (5-50 A) these concepts 
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are not well defined. Instead, in the molecular realm we distinguish between harmonic models, assuming atomic 
displacements follow Hooke’s law, and anharmonic cases, allowing for anharmonicities in the molecular force held. In 
harmonic models the current may show ballistic characteristics, with the conductance determined by the contacts only, 
or, “phonon tunneling”, when phonons off-resonance with molecular vibrations, cross the junction, showing features of 
quantum tunneling [^. We use Eq. (12) below to exemplify the concepts of ballistic and tunneling phonon dynamics 
in molecular heat conduction. Note that we use the notions of phonons and vibrations interchangeably, referring to 
the collective motion, excitations, of atoms. 


II. EXPERIMENTS: FOCUS ON HEAT CONDUCTION IN ALKANE CHAINS 

We begin with the presentation of several experiments, corresponding to the schemes in panels (a) and (b) of Fig. 
[l] This introduction serves to motivate modeling and computational works, but it does not aim to give readers a 
complete review over heat transfer measurements, provided e.g. in Refs. pm]. We focus on studies of length- 
dependent thermal transport in alkane chains with 2-25 methylene units. This molecule serves as a standard for 
phononic transport measurements for several reasons. First, these sigma-bond hydrocarbons are poor conductors of 
electricity, therefore thermal conduction takes place by phonons, not electrons. Second, the molecule is uniform and 
quasi one-dimensional. Lastly, at room temperature anharmonic effects are expected to be non-influential even in 
long chains of A ^ 100 units [35|. Altogether, alkane chains with N > 5 units are expected to conduct vibrational 
energy ballistically at room temperature, as confirmed by simulations |38j and experiments [inimMj- 

Solid-SAMs-solid junctions were studied in several works. The thermal conductance of Au-alkanedithiol-GaAs 
SAMs with 8-10 carbon units was measured in Ref. |31] with a 3a; technique, resulting in the room-temperature value 
(per unit area) of k ~ 25 — 28 MW m“^ K“^, independent of length. More recently, heat transport measurements 
with scanning thermal microscopy (SThM) were reported in Ref. [32j . with an alkane-thiol molecule adsorbed on a 
gold surface and a silicon tip approaching the SAM from above. The measured value for the thermal conductance is 
depicted in Fig. Sb). It shows a non-monotonic behavior: The conductance first rises with length for short chains, 
with the maximum conductance of 25 pW/K for a chain with four carbon atoms, but in longer chains k is reduced 
to ^ 10 pW/K. This crossover behavior appears in simulations, see Fig. [^a), to be discussed in more details in Sec. 

mm 

In ballistic transport the conductance reflects the molecule-contact interaction rather than intrinsic molecular 
properties. The effect of chemical bonding on interfacial heat transport in Au/alkane SAMs/Quartz junctions was 
examined in Ref. |33j . by studying chains with either methyl or thiol-termination, using the time-domain thermore¬ 
flectance (TDTR) technique, an optical pump-probe tool. The thiol-ending group forms a stronger (covalent) bond 
to the Au surface, relative to the weak (van der Waals) forces attaching a methyl group to a gold substrate. In 
accord with expectations, the interface thermal conductance was enhanced with the increase of bond strength: The 
interfacial thermal conductance (per unit area) measured for the thiol ending group was k=68 MW K“^. A 
lower value was reached for the interface missing thiol bonds, k = 36 MW m“^ K~^ [53]. In Ref. [35] it was also 
demonstrated experimentally that the thermal conductance of alkane SAMs can be systematically controlled-reduced 
by using metals with mismatched phonon spectra, e.g., Au and Pt. 

Bridged-type situations as depicted in Fig. [^b) were examined in several works, revealing information on mecha¬ 
nisms of heat flow from the dynamics of the transients. In one of the earliest studies, alkane chains were connected to 
large chromophores, azulene and anthracene, which effectively act as heat reservoirs [10]. The IVR process began by 
the electronic excitation of the azulene end-group, and after a fast internal conversion, vibrational energy localized 
on the azulene propagated through the bridge to the anthracene, with an IVR time constant of r ^ 5 — 6 ps. It was 
demonstrated that results were consistent with a ballistic transport model [lOj . At later times, equilibration with the 
surrounding solvent took place. 

A rather different approach, time-resolved flash-heating of SAMs on substrates, was utilized in Ref. [3D], see also 
[SHUT]. In this type of experiments, a thin metal film (gold) was heated by several hundreds of degrees within a very 
short time (^ Ips) using a laser pulse, to excite electrons near the metal surface. Subsequently, the lattice temperature 
rose due to electron-phonon coupling, and heat propagated from the base of the SAM, through the molecular chains, 
to the top (methyl) groups where it was detected with vibrational sum-frequency generation (SFG) spectroscopy. This 
signal is sensitive to disorder on the SAM, induced here by the temperature increase, as well as to frequency shifts. 
By varying the length of the alkane chain it was determined that heat flew in alkanes ballistically, with a velocity of 
9.5 A/ps, to yield a thermal conductance (per molecule) of 50 pW/K. 

Other recent studies further considered a bridged setup in solution and confirmed the ballistic energy transport 
mechanism at room temperature in alkane chains of 5-15 units, by using two-dimensional infrared spectroscopy (2DIR) 
[33][3D]. In this experiment, the energy transfer process was initiated by vibrationally exciting a “tag” group attached 
to one end of the molecule (e.g. N=N bond), while the vibrational energy arriving at the other end was recorded by 
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FIG. 2: Thermal conductance of alkane chains with increasing size, (a) Simulations under the harmonic 
approximation based on the Landauer formula, as reported in Ref. [38]. The reservoirs’ spectral functions are 
modeled by 7 ( 0 ;) = identical for L and i?, with Debye (cutoff) frequency u}c= 0.05 eV and a coupling 

energy a = 1 eV. Full line: T=50 K; dotted line: T = 300 K; dashed line: T=1000 K. (b) Measured values for the 
thermal conductance of SAMs of alkane thiols [HS-(CH 2 )Ar_i-CH 3 ] on gold substrates on Mica, using SThM with a 
silicon tip. Panel (b) is reprinted with permission from Reference |32j . Copyright 2014 American Physical Society. 


a “reporter” group, e.g., a carbonyl. Measurements at room temperature, in solution, revealed fast energy transport 
through alkanes with a speed of 14.4 A/ps and a mean free path of 14.OA, see Fig. Interestingly, while it was 
demonstrated that transport via the chain was ballistic, it was argued that heat propagated through the end-groups 
in a diffusive manner |35| . It was also demonstrated in Ref. [42] that by initiating dynamics with different tags, 
different chain bands (e.g., CC stretching, CH 2 twisting modes) contribute to the transfer process. 

Thermal transport experiments on hydrocarbons confirmed ballistic dynamics, the predominant mechanism in 
short-ordered systems. In contrast, diffusive dynamics is prevalent in complex systems (even of reduced dimensions), 
nanotubes and nanowires [smisj, polymers [S], peptides IHlISZj, as well as in small organic molecules |3S|. For 
example, in Ref. [37] diffusive transport of vibrational energy was observed in short 3io-helical peptides. The process 
was initiated by depositing excess energy e.g. onto a vibrational CD mode. Flow of energy through the helix was 
then detected optically by planting vibrational probes, isotopically-labeled CO modes, in several locations away from 
the source, serving as local thermometers. Measurements in Ref. [37] were consistent with a diffusion model, with a 
diffusivity constant of 2 A^/ps [Hj. 

In the following discussion over theoretical-computational frameworks we do not direct the issue of the ballistic- 
diffusive crossover in thermal conduction, as diffusive dynamics is not expected to fully develop in short molecular 
junctions, our focus here. More generally though we will evaluate the role of nonlinearities, deviations from the 
harmonic approximation for atomic displacements, in the heat transport behavior. What are the signatures of anhar- 
monicities in relatively small phononic conductors? First, one may observe the suppression of the thermal conductance 
with temperature, at high temperatures, when inelastic phonon scattering effects are at play |45| . Further, anhar- 
monicities can furnish nonlinear functionalities such as the thermal diode effect uni l46l l47] . which is specifically 
predicted to develop in the azulene-alkane-anthracene compound of Fig. [^b) at high temperatures [15] . Other dis¬ 
tinct signatures of vibrational anharmonicities are processes such as energy localization and unidirectional energy 
flow in molecules. Using IR Raman spectroscopy, it was manifested that in nitrobenzene vibrational energy could 
leave benzene-localized modes to occupy global modes, but energy could not propagate away after a nitro excitation 
[49H51| . Asymmetries in forward-backward transfer efficiencies certainly break the noninteracting phonon picture. 
We describe now quantum mechanical approaches for modelling vibrational energy flow in molecular junctions. 
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FIG. 3: Two-dimensional IR spectroscopy of vibrational energy trasnfer in alkanes, (a) Inferred energy transport 
time T^ax, from the tag mode to the reporter, and the tag-reporter cross-peak maximum amplitude Imax^ plotted as 
a function of chain length. The behavior of Tmax is approximated by a linear function, and from the (inverse) slope 
one receives the speed ^ 14.40 A/ps. The cross peak amplitude is approximated by an exponentially decaying 
function with a characteristic decay distance Rq ^ 14.6 A, representing the mean free path, (b) Schemes of energy 
transport mechanisms in alkane chains with different reporter groups (blue). Green arrows mark regions of ballistic 
propagation, red arrows mark (diffusive) steps from the tag group to the molecule, and from the molecule to the 
reporter. Yellow arrows point to the relaxation of the vibrational wavepacket. Also indicated are ballistic transport 
times, computed from the chain length divided by the energy transport speed 1440 m/sec. Reprinted with 
permission from Reference [3S]. Gopyright (2015), AIP Publishing LLC. 


III. METHODOLOGIES: HARMONIC APPROXIMATION AND THE LANDAUER FORMULA 

A. Formalism 

Under the harmonic approximation for atomic displacements, the molecule (“subsystem”) includes harmonically- 
connected masses, and it is linked linearly-harmonically to thermal environments composed of harmonic oscillators 
(“baths”). Since all terms are harmonic, one can diagonalize the model Hamiltonian to reach its fundamental normal 
modes (phonons) and solve the dynamics of this noninteracting problem exactly. The steady state transport behavior 
of the model can be arrived at by using different methods, including the generalized Langevin equation miSHlESHil] 
or, equivalently, the non-equilibrium Green’s function (NEGF) approach (211 UHl [SS] (derivations under the weak- 
coupling approximation can be found e.g., in Ref. [56].) Both methods yield a Landauer-type formula for the heat 
current, an analogue of the elastic-coherent expression for charge current m- 

Following Refs. [HIMIISSI, we outline the derivation of the Landauer heat transport expression from the principles 
of the Langevin equation approach. In this procedure, one begins with Heisenberg equations of motion for displacement 
and momentum operators, assuming a system-bath factorized initial condition at t = —oo with the reservoirs prepared 
in a thermodynamic equilibrium. It can be shown that after integrating out the baths’ coordinates, the subsystem’s 
displacements follow a generalized (quantum) Langevin equation in which the effect of the reservoirs is encapsulated 
within dissipation kernels and noise terms. These terms particularly depend on the baths’ spectral properties. Fourier 
transforming the generalized Langevin equation administers the steady state limit. In frequency domain, the linear 
equations for molecular coordinates can be readily solved to construct two-time correlation functions, to organize a 
closed-form expression for the heat current, a Landauer-type formula dHIMlIiS- For harmonic models, the quantum 
Langevin equation and the NEGF method can be further employed to yield closed expressions for the fluctuations of 
the heat current |58l |59| . 

We now provide the working expressions. The Hamiltonian for a harmonic junction reads 

H = Hs + Hl + Hr + Hi, (1) 

where the molecule Hr, the two contacts Hu, v = L,R, and the coupling term Hi are written as 

Ha = ^PaPa+^'>J'aKaUa, a = S,L,R 

^ ( 2 ) 

The vector Ua = gathers mass-normalized displacement operators for the three regions a = L, R, S. Similarly, 

Pa holds the conjugate momenta, T stands for matrix transpose. The force constants are organized into real symmetric 
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matrices Ks, and Kn, and the harmonic (separable) interaction of the molecule with the two reservoirs is given 
in terms of the force constant matrices K^s = [Asl]^ and Krs = Note that in general one does not need to 

assume a normal mode representation for the reservoirs as this may be achieved through diagonalization with unitary 
matrices Uu, 


u^K,u, = nl. 

In steady state, the phonon current across harmonic junctions is given by a Landauer-type expression m 

1 


3q = 


27r 


du) hioTioj) [nL{uj) — nii{ui)]. 


(3) 


(4) 


Here, n,^{uj) = — 1] ^ are Bose-Einstein distribution function for the v = L,R bath with the inverse temperature 

fiv = T{uj) stands for the transmission function calculated from [THl ICT] 


T{w) = Tr[G5(c.)ri(a;)Gg(a;)rfi(cc) 


(5) 


For small temperature differences, AT << T with AT = Tl — Tr, T = (Tr + Tr)/ 2 as the average temperature, the 
Landauer expression reduces to 


AT f°° 

jq = J dwhuTiio) 


dn(u}) 
dT ’ 


( 6 ) 


with n{u}) as the Bose-Einstein function at the averaged temperature T. The ingredients in Eq. ([^ are the Green’s 
function of the molecule, Go(a;) = (Gg(a;))^, and its self energy. 


GS(a.) = [(w + iQ+f -Ks- E£(a;) - EJ^(a;)] \ 

E:ic,) = [A,sfg:ic,)A,S. r,(..) = i[Ei;(a;)-E“(a;)], u = L,R 


(7) 


where r,y(w) is referred to as the spectral function for the leads. The unperturbed-isolated Green’s function for the 
baths are given by 


glit) = -e(t) gliu:) = H e^-^gl{t)dt = [(cc + * 0 +)^ - K,]-\ 

J—oo 


( 8 ) 


with 6{t) as the Heaviside step function. The Landauer equation @ relies on the harmonic approximation, yet even 
under this strong assumption it can be employed to address a broad range of problems: the behavior of heat conduction 
in one-dimensional chains, as opposed to the three-dimensional case, the role of disorder in thermal transport, as well 
as impurities, molecular structure (linear or T-shaped objects), reservoir’s spectral properties, and contact geometry. 

We now specify Eq. (H) and introduce a chain model with the first (last) atoms in the chain coupled to the left 
(right) baths. 


H = Hs + Y, 


IGL 


'i +I 


_ A;ni Y 

) 


+ E 

vGR 


Pr 1 2 


Ur — 


XrllN 


(9) 


The small-letter parameters. A, are elements in the matrices A, and respectively. The self energy matrices 
reduce to single elements, e.g., at the left contact 


[E£(a;)]„,„- = Eliuj)6n,n'SnA = - [ 

JO 




IGL 




( 10 ) 


The real part of the self energy presents frequency shifts to molecular modes, the result of the coupling to the bath. 
The imaginary part of the self energy is responsible for energy damping. Ignoring frequency shifts, it is convenient to 
express it in terms of a (real-valued) function 7 i/(w), defined from E0(a;) = It corresponds to the friction 

coefficient in the language of the quantum Langevin equation. 

Simulations of phononic heat current based on Eq. Q are feasible, but analytic results are limited to special cases, 
e.g, to the classical-high temperature limit of a uniform chain |19j . or to short atomic bridges |62] . To provide insights 
into the transport behavior, we exemplify Eq. 0 in the single-atom limit, taking into account a single oscillator 
which is directly coupled to the L and R reservoirs. We denote the frequency of this mode by ojo adhering to the 
notation employed in the literature, see e.g. Refs. [H31E11- The transmission function of this bridge is given by 






(uj^ - a;5)2 -h + 7it(w)]^ 


( 11 ) 
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FIG. 4: Quantum interference effects of phononic conduction in molecular junctions, (a)-(b) Examples of localized 
phonon modes (linear combinations of normal modes) in oligo(phenylene-ethynlene) (OPE3), a cross conjugated 
molecule. Numbering counts the eigenmodes. The modes may behave similarly to the Eano model (a) or the 
degenerate-levels model (b). (c) meta and para OPE3 configurations, (d) calculated transmission functions, and (e) 
the thermal conductance, considering out-of-plane modes only (z) or the complete spectrum involving (x, y, z) 
displacements. Reprinted with permission from Reference [70], Copyright (2013), AIP Publishing LLC. 


and the phonon current through the junction is 



diuHu} 


a;^7L(w)7i?,(w) 


(w 2 f -f a;2[7i(a;) -^ 7 ^( 0 ;)]' 


;[nL{uj) - UR^Uj)], 


( 12 ) 


see also Refs. jOH [OS] (albeit adopting a different notation). This simple expression describing quantum coherent 
transport allows us to clarify several concepts: (i) When incoming phonons are in resonance with the molecular mode, 
u! = ujQ, the transmission is only limited by contact effects. We refer to this type of motion as “ballistic” transport. 
Furthermore, in the ideal case with a perfect transmission, T{uj) = 1, we evaluate the conductance from ® and 
receive k = ^ dujHuj^ = k’^n'^T/ih, which is the (universal) quantum of thermal conductance IBTJl 1^ IB7] . 
confirmed experimentally in Ref. | 68 j . (ii) The damping function 71 ,, resulting from the molecule-bulk coupling, 
broadens the resonance allowing for “tunneling” of phonons, transmission of modes off-resonance with the molecular 
frequency, uj uiq. It can be shown that the tunneling contribution exponentially decays with bridge-length [38j . 
in analogy with the super-exchange mechanism of electron transfer [69] . (iii) The conductance grows with T at low 
temperatures, but once ksT > Hloq, it saturates as it approaches the classical limit. 


B. Applications 

The transmission function (§ can be evaluated using the so-called atomistic Green’s function approach [28] to 
provide the thermal conductance of a broad range of nanoscale objects, carbon nanotubes. Si nanowires, graphene 
sheet, graphite. In such calculations one typically computes the phonon dispersion based on a model for the interatomic 
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potential, then evaluates the transmission function employing a model for the self energies. Phonon transport in carbon 
nanotubes was examined in e.g. Refs. [niiis] using this analysis. 

We discuss now examples of Landauer-based calculations of heat transfer in relatively short molecules. Alkane chains 
and isotopically substituted disordered chains with 2-25 units were examined in Ref. [38) . showing the conductance 
to exhibit non monotonic features at room temperature with growing size, see Fig. This behavior results from 

the fact that the molecular vibrational spectrum may be largely altered with length for short molecules, as well as the 
degree of localization of some molecular normal modes. Since details of the molecule-heat baths (solids) coupling, the 
spectral properties of the solids, specifically the value of their Debye frequency, and the bulk temperature, determine 
which molecular modes participate in the conduction process, the conductance of ultra-small systems can be tuned 
with size, before ballistic motion takes over m- FigH a) displays results at different temperatures assuming ohmic 
thermal baths. Panel (b) depicts measurements of the thermal conductance in SAMs |32]j demonstrating a qualitative 
agreement in trends between theory and experiment. 

Reducing phononic heat conductance in molecular junctions (while maintaining high electronic conductance), is 
desirable for improving thermoelectric efficiency. Several design principles for molecules were tested in recent proposals, 
bringing modest success, as we discuss below. The challenge in the manipulation of phononic current is rooted in 
the bosonic statistics: Technically, if we replace the bosonic form in Eq. (§ by the corresponding fermionic function, 
the current would be determined by the transmission function only in the vicinity of the equilibrium Fermi energy. 
For molecular structures, the charge current then combines contributions from a single or few resonance levels. In 
contrast, the derivative of the Bose-Einstein distribution function in Eq. Q covers frequencies in a window up to the 
baths’ temperatures, potentially comprising substantial contributions from many molecular modes of different nature, 
localized and extended, to confuse dynamics. 

In Ref. [73] it was proposed that molecules made of two separate subunits with a weak chemical link, for example, 
TT-stacked aromatic rings, could show reduced thermal conductance relative to the case with a single unit, while 
maintaining good electrical conductivity (since electronic overlap through-space is preserved in tt—tt systems). Detailed 
simulations with first-principle parameters revealed a more intricate picture: The phonon conductance may indeed 
reduce (by about a factor of 2) in this design, relative to the case with a single unit, but in some cases, it increased 
due to an overall improvement of coupling in the junction. 

Reduction of thermal conductance in molecules by employing the effect of quantum interference was demonstrated 
in Ref. |70j . It was pointed out that in e.g. a benzene ring with a meta connection to the leads, the phononic 
transmission function (comprising in-plane vibrations) can exhibit destructive quantum interference features, resulting 
in the suppression of the thermal conductance compared to the linearly conjugated analogue, see Fig. This 
behavior is conceptually similar to the electronic case [74l[75]. However, while variations in electronic conductance 
due to interference effects may be of 1-3 orders of magnitude [741175], changes in phononic conductance due to quantum 
interference are rather modest, factor of 1.5 to 5 in Ref. m- As explained above, this is because phononic conductance 
collects contributions from many modes of different nature, whereas quantum interference effects influence individually 
specific modes. 


IV. METHODOLOGIES: ANHARMONIC INTERACTIONS 


The Landauer picture breaks down when deviations from the noninteracting normal-mode picture are significant. 
This is the case in systems with large-tunable anharmonicity [as in the FPU model m] and at high temperatures 
when inelastic phonon scatterings are influential. What tools are available beyond the Landauer-harmonic approach? 
Classical molecular dynamic simulations can be naturally applied beyond harmonic force fields, to include interactions 
to all orders [a EH Ell [78]. Molecular dynamics simulations of heat conductance in simplified models revealed the 
role of disorder and reduced dimensionality |7MST] as well as anharmonicities on transport mechanisms and 

the occurrence of nonlinear phenomena such as negative differential thermal conductance [87j . thermal rectification 
[HHlIHni, gating EO], and memory elements El- 

Quantum mechanical treatments are typically limited to a certain range in parameters. Among relevant methods we 
mention the non-equilibrium Green’s function technique, which is perturbative in the nonlinear interaction strength 
[27l[28], and master equation approaches, which may include the effect of anharmonicity exactly (as long as the model 
is minimal), but are often limited to models with weak system-bath couplings |45l|63l|92|. Other schemes are based on 
time-scale separation between slow and fast modes |S3|, phenomenological self-consistent treatments [iniiii[Miii5|, 
and numerically-exact methodologies, computationally limited to minimal modes |96H98j . 
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(a) 





FIG. 5: (a) Illustration of the SCR model in a chain model with N = 5 atoms. Under the self consistent conditions, 
JL = jn- The conditions j 2 = ja = jd = 0 set the temperatures Ti, i = 2 — 4. (b) An example of a spatially 
asymmetric model (the different sizes reflect different masses), which realizes a thermal diode effect - only in the 
quantum regime [94]. (c) Graphene nano constriction and (d) its self-consistent temperature profile (K) in the 
quantum SCR model. Temperatures were determined self-consistently for the gray atoms in the central region. A 
classical approximation provided very similar results. Panels (c) and (d) are adapted with permission from Ref. [95j . 

Copyrighted 2013 American Physical Society. 


A. Self-Consistent Reservoirs 

The self-consistent reservoir (SCR) method introduces effectively anharmonic interactions into harmonic models by 
allowing the system to exchange energy with internal reservoirs which mimic many-body interactions. In Fig. [^a)-(b) 
we exemplify the model with a one-dimensional chain. The system includes N atoms (not necessarily identical) in a 
junction setting, connected to each other by harmonic springs. In addition, each particle is coupled to an independent 
harmonic (Langevin) thermal reservoir, also referred to as a “thermal Biittiker probe” [99] . The temperatures at 
the “exterior” reservoirs, attached to the first and last particles, are set at and Tr, but the temperatures of the 
internal baths are determined in a self consistent manner, by demanding that in steady state, on average, there is no 
net heat current from these baths to the physical system. 

The SCR model has been investigated intensively with toy models. It was originally proposed in the classical domain 
[1001 llOlj , to prove the emergence of diffusional dynamics and Fourier’s law of heat conduction upon the application of 
the interior reservoirs [102] . The model was later extended to describe heat conduction in quantum systems, in linear 
response [1^110311104] and beyond that, analytically-approximately [10511106] . and numerically-exactly [94l|95]. What 
are the signatures of anharmonicity within the SCR model? The interior reservoirs can gradually tune the transport 
dynamics, from ballistic-elastic to diffusive. Further, the quantum SCR method can support a diode behavior (which 
is identically missing under the Landauer expression, as well as in the classical SCR method), due to the interplay of 
quantum effects, spatial asymmetry, and effective anharmonicity, see a possible setup in Fig. [^b) following Ref. [94] . 

Beyond toy models, the SCR method has been recently applied in Ref. [95] for the study of quantum thermal 
transport in nanostructures such as the two-dimensional graphene constriction depicted in Fig. [^c). The temperature 
profile in the central region was evaluated with a quantum SCR treatment, displayed in Fig.lmd). Interestingly, it 
was found that the quantum profile closely matches classical results, suggesting that low-frequency modes, for which 
quantum and classical statistics agree, largely contribute to transport in this system. 

In conclusion, while missing genuine anharmonic interactions, the SCR method offers a simple mean for study¬ 
ing phononic conduction in nanostructures while including both quantum effects and (effective) inelastic phonon 
scatterings. 


B. Non-equilibrium Green’s Function Approach 

1. Formalism 

The non-equilibrium Green’s function technique offers an elegant-powerful mean for computing equilibrium and out- 
of-equilibrium properties of many-body systems [MHiini. For quantum transport problems, the NEGF approach 
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provides a computational framework for incorporating interactions (anharmonicities, in the heat transport problem), 
going beyond the Landauer expression. We begin by introducing the Meir-Wingreen formula for the steady state heat 
current flowing through an anharmonic system which is linearly (harmonically) coupled to two harmonic reservoirs L 
and j?. 


3q = 


1 

— da;fiwTr[G<(a;)E>(w)-G>(w)E<(a;)]. 

27r Jq 


(13) 


Here, 'E^{ui) = —inL,{uj)T^(uj) and Ej)’(a;) = —i[l + n^(uj)]T^{u!), are the self-energy components for u = L,R. 
The harmonic part of the Hamiltonian was described in Sec. |HI A[ Adding nonlinear terms to Hs, the molecular 
Green’s function G now incorporates anharmonicities. One can formally write down a closed set of equations for its 
components in terms of the nonlinear self-energy E„, 


G’- = [(G5)-i - EO] ^ = [(a; + zO+f 
G< > = (E<’> -H E<’> -f E<’>) G“. 


-1 


(14) 


For simplicity, we suppressed the w dependenc e of these functions. Go in Eq. (14) precisely corresponds to the Green’s 
function defined in the harmonic case in Sec. |HI[ The second equation in (14) is known as the “Keldysh equation” 
[109) . In the absence of anharmonic interaction, E„ = 0, the Meir-Wingreen expression reduces to the Landauer 
formula with the transmission function Eq. (§. In the interacting case it can be organized and expressed in terms of 
an effective, temperature-dependent, transmission function 1113) . 

Tes{uj,TL,TR) = Tr[G”rLG“rfl] 

^ -Tr[G'T„G“(niri - nflrfl)-iG'-E<G“(ri - Tr)], 


+ 


2{nL-nR) 


(15) 


where F^ = i(E” — E“) and rii, denotes the Bose-Einstein distribution function of the v bath. Under a symmetric 
coupling, Fl = Ffl = F, the transmission function reduces to 


T:r{u;,TR,TR) = iTr[A(a;)r(cc)], 


(16) 


with A(w) = i[G’'(w) — G“(cli)] as the spectral function of the molecule. 

The Meir-Wingreen formula (131 is formally exact, but approximations should be employed to compute the nonlinear 


self-energy, hence the current. Computational treatments include the standard quantum field theory-type perturbative 
scheme using the Feynman diagrammatic approach [113H115] . and self-consistent approaches ITTBIIl 19) in which 

the harmonic Green’s function Go in E„ is replaced by the full nonlinear G. In fact, the self-consistent approach 
under certain conditions provides results in agreement with quantum master equation treatments [12011121) , when 
the molecule-reservoir coupling is sufficiently weak. 

What is the physical content of the retarded nonlinear self-energy E”? Its real part is responsible for frequency 
shifts of molecular vibrational modes. The imaginary part corresponds to the finite life-time of phonons, eventually 
responsible for the development of diffusive dynamics. It should be mentioned though that to fully realize the diffusive 
regime, high order phonon scattering processes should be included in the nonlinear self-energy. Once extracting the 
phonon lifetime Tq, the thermal conductivity (denoted here by k) can be computed from the Boltzmann transport 
equation or more simply, using the standard kinetic theory formula k = where Cq is the heat capacity 

per unit volume of mode q and Vq is the phonon group velocity [122) . 


2. Application 


The Landauer formula @ was derived under the harmonic approximation thus it rules out thermal rectification 
effects and other nontrivial nonlinear functionalities m- The NEGF formalism, in contrast, allows us to evaluate the 
role of a nharmonicities in energy transport. We exemplify the utility of the method on the single atom junction [recall 
Eq. (^llb], allowing the oscillator to further interact with a quartic pinning potential with the Hamiltonian = f^o 
[123] .^o the lowest (first) order in e, the components of the no nlin ear self-energy are given by E”(a;) = E“(a;) = 


3zIieGQ (t = 0), E<(a;) = E>(a;) = 0. Using these terms in Eq. (15), the transmission function with nonlinearities 


Tn{to) can be expressed in terms of the harmonic result T{ui) [124) . 

Tn{ta-,TL,TR) = {1 + Aiuj-,TL,TR))T{uj), 

A{u;-,Tl,Tr) = 6iheG<{t = 0)Re[G5(w)]. 


( 17 ) 
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FIG. 6: Self-consistent NEGF calculations of the thermal conductance as a function of temperature T for an 
anharmonic junction. The single harmonic oscillator junction (e = 0) is supplemented by a quartic onsite potential 
Hn = fuQ. e is given in units of eV/ (amu^A^). The reservoirs are defined by semi-infinite chains of harmonic 
oscillators with the nearest-neighbors force constant k = leV/ (amuA^) and an onsite element ko = 0.1/c. The 
“molecule” includes a single oscillator with Kg = 1.1k, and it is coupled to the left and right solids through 
A:^;f Q = Aq^ = —0.25/c. The inset compares the thermal conductance as calculated using the self-consistent 
procedure (solid line) and the first-order expression ( [l7| (dotted line) for e = 2eV/ (amu^A^). Adapted with 

permission from mi¬ 


lt should be recognized that the lesser component of the Green’s function, G(5'(t = 0), depends on temperature. As 
a result, the effective transmission function becomes temperature dependent and the junction can manifest nonlinear 
functionalities such as the diode effect. 

Fig. (§ displays the thermal conductance of this single anharmonic oscillator junction, plotted against the averaged 
temperature T. Simulations were performed by employing the self-consistent NEGF approach, replacing Gq by the full 
G in E„, and solving the Keldysh equation (141 iteratively, to compute the current from Eq. (13) |123j . It is evident 
that at low temperatures the junction’s thermal conductance is dominated by elastic-harmonic forces. In contrast, at 
high temperatures the thermal conductance is significantly reduced when the anharmonic pinning potential is turned 
on, as inelastic phonon scattering processes largely influence transport. 

The NEGF technique as described here is limited to small-atomistic models. Efforts are dedicated to link and 
interface NEGF results with coarse-gained modeling such as the Boltzmann transport equation. This would allow 
simulations of larger nanoscale systems, and the resolution of e.g. ballistic-diffusive crossover with parameter-free 
quantum mechanical first principle approaches m- 


C. Quantum Master Equation 

Projection Operator approaches, such as the Nakajima-Zwanzig projection operator technique [125j . trace out 
the dynamics of the (non-interesting) bath degrees of freedom, to reach equations of motion for relevant-subsystem 
coordinates, typically, under the assumption of fast dynamics for the bath relative to the subsystem. At the heart of 
such approaches is the conceptual separation of the model into three constituents: subsystem of interest, environment, 
and an interaction term between the components. Quantum master equations (QME) which describe the reduced 
dynamics of the density operator are immensely useful for modelling electron, proton, and exciton transfer, vibrational 
relaxation and spin dynamics in condensed phases |126j . As well, in the related optical QME the radiation field plays 
the role of a thermal environment [125) . Standard approximations involved in the framework of QME are: (i) weak 
subsystem-bath coupling, allowing a perturbative expansion in the interaction strength, (ii) Markov approximation. 
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Spin-boson model Extended models 


FIG. 7: Thermal nanojunctions in the energy representation for the molecular subsystem, (a) Generic anharmonic 
junction, (b) Damped harmonic oscillator model, generalized to the non-equilibrium case, (c) Spin-boson model 
with an anharmonic, two-level system, coupled to harmonic baths, (d) An extended-site model which can describe 
e.g. spin and exciton transfer in chains. The arrows describe an energy transfer process, with a certain site being 

excited simultaneously with the de-excitation of a different site 


assuming short-lived memory effects in the bath, and (iii) rotating wave approximation (often employed in the quantum 
optics literature), leaving out off-resonant terms. 

QME traditionally communicate the dynamics, population and coherences, of a subsystem of interest. More recently, 
efforts were put forward for the development of related schemes for transport coefficients, to calculate (charge, heat, 
spin) currents in non-equilibrium situations. Advantages of QME over other approaches such as the NEGF are: (i) 
It can treat exactly intrinsic anharmonicities, yet obviously, limited to small “impurity” models, (ii) Analytic results 
can be derived, to expose the role of different factors in the heat transport behavior. The QME approach can be 
conveniently formulated in the energy basis of the subsystem. Some examples in the context of vibrational heat flow 
in molecular junctions are illustrated in Fig. 


1. Weak molecule-reservoir coupling: The Bloch-Redfield-Markov equation 

We bring here a brief description of QME for vibrational heat transport in junctions, under the assumption of weak 
molecule-bath coupling, by following Refs. [531112]- Fig- [^a) displays the model, an Wstate subsystem (possibly 
many-body nonlocal states) and two heat baths. Excitation and relaxation processes in the subsystem are induced 
by complementing processes in the reservoirs. The total Hamiltonian is 

H = Hs + Hl + Hr + Vl + Fr, (18) 

where 

Hs = Y,En\n){n\, = v = L,R. (19) 

n 

The energy spectrum of the subsystem (molecule) may be derived from a truly anharmonic potential without 
further approximations. In the damped harmonic oscillator model of Fig. En = nhuiQ. The molecular system is 

coupled to thermal baths u = L,R, not necessarily harmonic [see examples in Refs. [B31|33]], with the subsystem 
operators S" = J2m n the operator B", given in terms of the v bath degrees of freedom. Leaving out 

details, it can be shown that the population of the subsystem eigenstates obey kinetic equations [31] 

Pn{t) = -Pn{t) Y. \SlSK^m + E IIW , (20) 

with rate constants 



k: 


( 21 ) 
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This result is correct to the lowest nontrivial order in the subsystem-bath coupling. It was derived under a Markovian 
approximation, and after decoupling the diagonal and off-diagonal elements of the reduced density matrix. The 
operators in Eq. are given in the interaction representation, (...)t represents an average with respect to the 

equilibrium distribution of the two baths, Em,n = Em — En- Applying similar steps, the steady state heat current 
across the system [depicted in Fig. [^a)] explicitly follows [531 [S2] 


Jq = 


- y 

n'ym 


Em,n\St,n?[Pnki 


-Pr 


-\y Em,n\Sy?[Pnk^_ 


— p y 


( 22 ) 


Here, the population are the long-time solution of Eq. (20 1 . This symmetric expression is quite intuitive: The 
current is given by calculating the net transfer process from all transitions between the n and m states. Rate 
constants are multiplied by energy differences and weighted by the population of states. A more general result, 
involving contributions from coherences, was derived in Ref. [45j . 

We exemplify this expression in the “atomic limit”, with a junction comprising of a single oscillator. In the 
harmonic limit, using the bosonic creation and annihilation operators, ^ Hs = fnooblbo, 






{bi g + bv^q){b\ -I- 6 o), one can derive the Landauer-like result, 


jq = hWQ 


7L(^o)7fl(^o) 

7l(wo) + 7fi(wo) 


[nL(a;o) - nR{oJo)] , 


(23) 


where 'ji^{uj) = t: 12 J2q ^(5(a; — ujq) stands for the mode-bath coupling energy, consistent with the definitions of the 
damping function in Sec. Ill A Eq. (23) can be obtained from Eq. (12) taking the limit 'ji, wq, and it describes a 
resonant-sequential transmission process, as only reservoirs’ modes which match the molecular frequency ujq scatter 
through the junction. 

Beyond the harmonic approximation, the “spin boson model” is a simple yet truly nontrivial example for anharmonic 
heat conduction. In this case, the subsystem spectrum is truncated to include only two states with a gap ujq, see Fig. 


7 c). 


The spin-boson Hamiltonian includes a two-state “impurity”, Es = interacting with harmonic baths 

Hi, = hiUqbl gb,^^q according to 14 = i^l q + ^v,q)- Here Ux and Uz are the Pauli matrices. Using Eqs. 

([2^- ’ ’ 


I ',o 

one receives now the heat current 0 


jq = bwi 


1l{i^o)ir{i^o) 


7l(wo)[1 -I- 2nL{uJo)] + 7fl(wo)[l -I- 2ni^(a;o)] 


[nL(a;o) - n_R(wo)]. 


(24) 


Since the baths temperatures enter the transmission probability here, this junction can display the thermal diode 
effect when spatial asymmetries are presented [46j . This nonlinear functionality is a definite fingerprint of anharmonic 
terms. 

The weak-coupling second order QME presented here can be extended systematically to higher orders, to account 
e.g. for co-tunneling processes m- Insight can be also earned by writing down Fermi-golden rule rates for multi¬ 
phonon scattering processes, then constructing the overall heat transfer rate, or the thermal current, based on the 
framework of kinetic master equations |S1 112811129] . For a cubic anharmonicity, for example, Fermi golden rule rates 
consist two types of (energy conserving) processes: a vibrational mode may decay into two others, or two modes can 
collide to produce a third mode. This type of analysis is particularly useful for simulating vibrational energy transfer 
in amorphous materials, glasses m, proteins |129L ll31H133j and across interfaces |3H] • 


2. Strong molecule-surface coupling: polaronic equations 


Equation (22) was derived under the assumption of a molecule weakly attached to the thermal reservoirs. It 
predicts a linear enhancement of the thermal conductance with increasing coupling energy 7 , see for example Eq. 
(23). Obviously, at strong coupling, corresponding to a molecule tightly attached to the interface, this trend should 
break down. By using the so-called non-interacting-blip approximation (NIBA) |134| . which is related to a polaronic 
treatment mi: it was recently demonstrated that the thermal conductance of the anharmonic spin-boson model 
(with an ohmic dissipation) displays a crossover behavior with the system-bath coupling parameter, satisfying k oc 

* ... / . 

] a = ul -\- api ISS]. We used the ohmic damping function, 71 ,(w) = Tra,yUJe with ujc the Debye 


7 2 / ,- \ 1 —2q; 

/ hujc \ 


yksT I 


(cutoff) frequency of the phononic baths and a a dimensionless parameter. It is also useful to define Er = 2a,^ojc 
as the reorganization-interaction energy of the molecular mode with the baths’ phonons. This NIBA result for the 
conductance agrees with numerically exact Monte-Carlo simulations, at high temperatures EZ]. 
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FIG. 8: Heat transfer in an anharmonic nanojunction, the spin-boson model of Fig. [^c). (a) Current as a function 
of the molecule-bath coupling parameter a = + C(r, demonstrating the breakdown of the linear scaling predicted 

by the weak coupling Bloch-Redfield-Markov QME [35], with a suppression of the current at strong coupling. Panel 
(b) Zooms over the weak coupling regime where different approximations agree. = 2a;o and ksTpi = cjq, 

Wc/wq = 10. (c) Thermal conductance k = jq/AT as a function of temperature in the linear response regime, 

AT = Tr — Tfj. NIBA (polaronic QME) correctly behaves only at high temperatures. Wc/wq = 10 and 

CiL.R = 0.0172. For more details, see Ref. [35] . 


Fig. [^displays the heat current in the spin-boson model, a representative of an anharmonic junction, under different 
approximation schemes: a standard weak-coupling treatment (Bloch-Redfield-Markov QME), NIBA-polaronic QME, 
an NEGF-Redfield approach |112 | and a more accurate variant of NEGF m. as well as numerically-exact influence 
functional path integral simulations |138j , performed based on the mapping of the model onto a fermionic description. 
For more details, see Ref. |35|. We make the following observations: (i) The trend jq cx a fails beyond the very- 
weak coupling limit, (ii) At high temperatures, ksT > Hujq, the thermal conductance decays with the increase of 
temperature, a manifestation of the intrinsic anharmonicity of the junction. This trend is correctly captured by QME 
calculations, see Fig. [^c). In contrast, at low temperatures this trend is reversed since thermal occupation factors of 
bath modes dominate the behavior, rather than temperature-dependent (anharmonic) effects on the junction. NIBA 
simulations, however, fail to capture this turnover behavior. Approximate techniques which can bridge between 
different regimes, to consistently describe transport e.g. at weak and strong coupling, are of great interest m , as 
well as the derivation of bounds for heat transport in anharmonic junctions |140j . 


D. Numerically-Exact Simulations 

The methods detailed so far are limited to describe heat conduction in certain physical regimes due to their underly¬ 
ing approximations, as exemplified in Fig. [^ It is highly desirable thus to develop first-principle techniques which can 
go beyond perturbation theories, to reveal the crossover between different regimes and to provide benchmarks for ap¬ 
proximate theories. Obviously, such nonperturbative treatments are numerical in nature and computationally limited 
to the simulation of minimal nanojunctions such as the two-bath spin-boson model. The dissipative dynamics of a two- 
state system coupled to a single heat bath, the (original) spin-boson model, has been examined in many studies |141j . 
Among the numerically exact machinery employed to treat it we list the quasi-adiabatic propagator path-integral 
(QUAPI) approach [14211143] . which builds upon the Feynman-Vernon influence functional, the multiconfiguration 
time-dependent Hartree (MCTDH) method, a wavefunction theory |144j . and Monte Garlo (MG) simulations |145) . 
These tools were developed to hand over the dynamics of the reduced density matrix of the subsystem. However, 
a current operator is a more complex object in nature since it is explicitly defined in terms of the baths degrees of 
freedom. 

Generalizations of numerically exact tools to simulate phononic heat current were reported in several recent studies: 
The multilayer multiconfiguration time-dependent Hartree theory was extended in Ref. [96| , and it revealed a turnover 
behavior of the current with increasing coupling strength, similarly to Fig. [^a). Influence-functional path integral 
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simulations were employed in Refs. |98l I146j . by mapping boson environments into fermion baths |138j . Monte Carlo 
simulations of the (linear response) thermal conductance were performed in [97) . by mapping the spin-boson model 
onto an Ising model with long-range exchange interactions. More recently, influence function path integral expressions 
were generalized to describe heat exchange between a subsystem and a heat bath |147) . though a calculation of the 
steady state heat current under a QUAPI-type approach is still missing. 

These numerically exact treatments are limited to a minimal anharmonic heat conducting nanojunction, a qubit, 
yet they are immensely useful for our understanding of quantum transport in hybrid junctions m in a broad range 
of parameters: under weak/strong coupling of the subsystem to the heat baths, low/high temperatures, slow/fast 
subsystem dynamics relative to the baths’ motion, and when different baths are employed, e.g., with rich spectral 
properties, possibly taking into account bath anharmonicities. Future studies will aim to bridge the gap between 
minimal models and more physical realizations of anharmonic junctions. 


V. CONCLUSIONS 

In this review, we focused on the problem of vibrational energy transfer (or thermal/heat/phononic conduction) in 
a junction geometry with a molecule bridging two heat reservoirs. We described recent experiments on vibrational 
heat transfer in molecular junctions focusing on alkane chains. We then reviewed methodologies useful for calculating 
thermal conduction in the quantum domain, starting from coherent-harmonic phononic conduction in nanojunctions, 
and concluding with highly-anharmonic yet minimal qubit constructions. We began by recounting the Landauer ap¬ 
proach, appropriate for the study of elastic transport. This simple tool can be readily employed to simulate phononic 
conduction in extended systems, as long as the molecular normals modes and their hybridization with the contacts are 
given. The Landauer-elastic scattering picture breaks down when interactions (anharmonicities) are presented. We 
outlined several techniques which account for anharmonic effects: The self-consistent reservoirs method mimics inelas¬ 
tic scattering of phonons with temperature probes. The non-equilibrium Green’s function technique as described here 
accounts for genuine many-body effects, albeit perturbatively, and it is limited to simple-minimal models. Reduced 
density matrix projection operator techniques can be extended to study thermal conduction of “impurity models”, 
when a molecule is strongly coupled to macroscopic solids. Finally, numerically exact techniques can simulate heat 
exchange in minimal models such as the spin-boson nanojunction, providing benchmarks for approximate techniques. 

Developing techniques for computing quantum heat transfer in nanostructures, bridging the gap between minimal- 
model junctions and real systems, is highly desirable. Ongoing efforts are dedicated towards the development of 
semiclassical techniques in which quantum statistics is implemented within classical molecular dynamics simulations 
[14811149| , construction of quantum corrections to classical expressions [251 1150j , and computations using Boltzmann 
transport equation with input from first-principle atomistic approaches [122j . 

In this review we only considered conduction of heat due to vibrations of nuclear degrees of freedom. However, 
when a molecule links metal electrodes, electrons as well considerably participate in the heat conduction process. 
Understanding how the interaction between electrons and atomic vibrations affect electron dynamics as well as energy 
transfer in real systems is a fundamental-formidable task, groundwork for making further progress in molecular 
nanotechnology. 
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